<h1 id="线性方程的迭代方法">线性方程的迭代方法</h1>
<p>本文将介绍求解线性方程 $Ax = b$ 的几种方法。 其中 A 是一个大型稀疏矩阵，通常以算子的形式给出，例如在偏微分方程中。这类问题的规模太大，以至于像LU分解之类的直接方法受内存限制无法使用，故需要使用迭代求解。</p>
<p>[toc]</p>
<h2 id="静态stationary方法">静态（Stationary）方法</h2>
<p>采用如下迭代格式： $$ x_{k+1}=-M<sup>{-1}(A-M)x_k+M</sup>{-1}b\tag{1} $$ M是一个相对于A更便于求逆的矩阵。可以很方便的验证上述格式在$Ax=b$ 是不动点。M有3种常用的选择</p>
<ol>
<li>$A$ 的对角元 $D$ (Jacobi 格式)，其具体迭代为： $$ x^{(k+1)}<em>i={1\over a</em>{ii}}\left(-\sum_{j\ne i}a_{ij}x^{(k)}_j+b_i\right)\tag{2} $$</li>
<li>$A$ 的下三角元 $D+L$ (Gauss-Siedel 格式)，其具体迭代为： $$ x<sup>{(k+1)}<em>i={1\over a</em>{ii}}\left(-\sum_{j&lt;i}a_{ij}x</sup>{(k+1)}<em>j- \sum</em>{j&gt;i}a_{ij}x^{(k)}_j+b_i\right)\tag{3} $$</li>
<li>以上二者线性组合 ${1\over \omega}D+L$ (SOR 格式)，其具体迭代为： $$ x_i^{(k+1)}={\omega\over a_{ii}}\left(b_i-{\omega-1\over\omega}x_i<sup>{(k)}-\sum_{j&lt;i}a_{ij}x_j</sup>{(k+1)}-\sum_{j&gt;i}a_{ij}x_j^{(k)}\right)\tag{4} $$</li>
</ol>
<h2 id="共轭梯度法conjugate-gradient-cg">共轭梯度法(Conjugate gradient, CG)</h2>
<p>当$A$ 对称正定时，可将原问题化为最小化 $\phi(x)={1\over 2}x<sup>TAx-b</sup>Tx$. 从一个点 $x$开始，沿一个方向 $p$ 搜索 $\phi(x+ap)$ 的最小值，得到 $$ a=\arg\min\phi(x+ap)={p^Tr \over p^TAp }, \quad\text{where }r=b-Ax\tag{5} $$</p>
<p>再把 $x+ap$ 当做新的起点，换个方向重复以上步骤就可以不断减小 $r$ 直到找到解。剩下就是选择搜索的方向。注意到$\nabla\phi(x)=-r$, 即 $r$ 是 $\phi(x)$ 下降最快的方向，故可以选取 $r$ 作为每次搜索的方向。但实际上有更好的办法。记 $P_{i-1}=[p_0, \cdots, p_{i-1}]\in\mathbb{R}^{n\times i}$ 为前 $i$ 次的搜索方向，则 $x_i=x_0+P_{i-1}y_{i-1}$, 其中 $y_{i-1}$ 为每次迭代的 $a$ 。则在第 $i$ 次迭代中 $$ \begin{aligned} \phi(x_i+ap_i)&amp;=\phi(x_0+P_{i-1}y_{i-1}+ap_i)\ &amp;=\phi(x_0+P_{i-1}y_{i-1})+ap_i<sup>TA(x_0+P_{i-1}y_{i-1})+{1\over2}a</sup>2p_i<sup>TAp_i-ab</sup>Tp_i\ &amp;=\phi(x_0+P_{i-1}y_{i-1})-ap_i<sup>Tr_0+ay_{i-1}</sup>TP_{i-1}<sup>TAp_i+{a</sup>2\over2}p_i^TAp_i \tag{6} \end{aligned} $$ 如果 $P_{i-1}^TAp_i=0$ , 即 $\forall j&lt;i,p_j^TAp_i=0$ ，则上式变为 $$ \phi(x_i+ap_i)=\phi(x_i)-ap_i<sup>Tr_0+{a</sup>2\over2}p_i^TAp_i\tag{7} $$ 注意到现在 $a$ 和 $y_{i-1}$ 已经分开了。则现在有 $$ \underset{y\in\mathbb{R}<sup>{i+1}}{\arg\min}\space\phi(x_0+P_iy)=[y_{i-1},a]</sup>T \tag{8} $$ 设 $\mathcal{K}<em>i=\left&lt;p_0,\cdots,p_i\right&gt;$ ,其中尖括号表示向量张成的空间，$\mathcal{K}<em>i+x_0={x_0+x|x\in\mathcal{K}<em>i}$ 。则按这种方法迭代的 $x_i$ 满足 $$ x_i=\underset{x\in x_0+\mathcal{K}</em>{i-1}}{\arg\min}\space\phi(x)\tag{9} $$ 由此，我们选择 $p_i$ 为 $r_i$ 减去其在 $A\mathcal{K}</em>{i-1}={Ax|x\in\mathcal{K}</em>{i-1}}$ 上的投影。记 $(A\mathcal{K}<em>{i-1})^\perp={x|\forall x'\in\mathcal{K}</em>{i-1},x^TAx'=0}$, 则有 $$ r_i=\gamma_ip_i+z_i,\quad\text{where }z_i\in A\mathcal{K}<em>{i-1},\space p_i\in(A\mathcal{K}</em>{i-1})^\perp\tag{10} $$ 到此我们已经得到了一个算法。但其在每步中需要使用 Gram-Schmidt 方法来求 $p_i$ 。下面说明 $p_i$ 是 $r_i$ 和 $p_{i-1}$ 的线性组合，从而对算法优化。先证几个命题。</p>
<p><strong>Proposition 1:</strong> $\mathcal{K}_k=\left&lt;r_0,Ar_0,\cdots,A^kr_0\right&gt;$</p>
<p><strong>Proof:</strong> 使用归纳法。$k=0$ 时显然成立，设命题对于 $k$ 成立，注意到 $$ r_i-r_0=-AP_{i-1}y_{i-1}\in A\mathcal{K}<em>{i-1}=\left&lt;Ar_0,\cdots,A^ir_0\right&gt;\tag{11} $$ 又 $z_i\in A\mathcal{K}</em>{i-1}$, 则 $$ \begin{aligned} p_i={1\over\gamma_i}[r_0-(z_i+r_0-r_i)]\in\left&lt;r_0,Ar_0,\cdots,A<sup>ir_0\right&gt;\ \mathcal{K}<em>i=\left&lt;\mathcal{K}</em>{i-1},p_i\right&gt;=\left&lt;r_0,\cdots,A</sup>ir_0\right&gt;\quad\quad\quad\quad\square\tag{12} \end{aligned} $$ <strong>Proposition 2:</strong> $P_{i-1}^Tr_i=0$</p>
<p><strong>Proof:</strong> 由前得知，$x_i=x_0+P_{i-1}y_{i-1}$ , $$ \begin{aligned} y_{i-1}&amp;=\arg\min\phi(x_i)\ &amp;=\arg\min\left{\phi(x_0)+x_0<sup>TAP_{i-1}y_{i-1}+{1\over2}y_{i-1}</sup>T P_{i-1}<sup>TAP_{i-1}y_{i-1}-b</sup>TP_{i-1}y_{i-1}\right}\ &amp;=(P_{i-1}<sup>TAP_{i-1})</sup>{-1}P_{i-1}<sup>Tr_0\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad(13)\ P_{i-1}</sup>Tr_i&amp;=P_{i-1}^T(r_0-AP_{i-1}y_{i-1})=0\qquad\qquad\square \end{aligned} $$</p>
<p><strong>Proposition 3:</strong> $\forall j\ne i,\space r_i^Tr_j=0$</p>
<p><strong>Proof:</strong> 由命题二得 $\forall x\in\mathcal{K}<em>{i-1},\space x^Tr_i=0$ ，故只需证明 $r_i\in \mathcal{K}<em>i$ 。使用归纳法，当 $i=1$ 时显然，若对 $i-1$ 成立，则 $r_i=r</em>{i-1}-a</em>{i-1}Ap_{i-1}$ ，又由命题一得 $A\mathcal{K}_{i-1}\sub\mathcal{K}<em>i$，即 $Ap</em>{i-1}\in\mathcal{K}_i$，则 $r_i\in\mathcal{K}_i\qquad\square$</p>
<p>由命题三可知，$\mathcal{K}_i=\left&lt;r_0,\cdots,r_i\right&gt;$, 且 $r_0,\cdots,r_i$ 构成一组正交基。</p>
<p><strong>Proposition 4:</strong> $p_i\in\left&lt;r_i,p_{i-1}\right&gt;$</p>
<p><strong>Proof:</strong> 首先，$z_i$ 满足 $$ z_i=\underset{z\in A\mathcal{K}<em>{i-1}}{\arg\min}\space\left|r_i-z\right|^2\tag{14} $$ 可设 $z_i=AP</em>{i-1}y_{i-1},\space y_{i-1}\in\mathbb{R}^i$，则 $$ y_{i-1}=\arg\min \left|r_i-AP_{i-1}y_{i-1}\right|^2\tag{15} $$ 设 $y_{i-1}=[w_{i-1},\space\beta_{i-1}]^T$，又 $r_{i-1}-r_i=Ap_{i-1}$，易得 $$ \begin{aligned} \left|r_i-z_i\right|^2&amp;=\left|\left(1+{\beta_{i-1}\over a_{i-1}}\right)r_i-{\beta_{i-1}\over a_{i-1}}r_{i-1}-AP_{i-2}w\right|^2\ &amp;=\left(1+{\beta_{i-1}\over a_{i-1}}\right)<sup>2\left|r_i\right|</sup>2+ \left({\beta_{i-1}\over a_{i-1}}\right)<sup>2 \left|r_{i-1}+AP_{i-2}\left({a_{i-1}\over\beta_{i-1}}w_{i-1}\right)\right|</sup>2\space(16) \end{aligned} $$ 最后一式第二项同 $\left|r_{i-1}-z_{i-1}\right|^2$ 比较得 $$ w_{i-1}=-{\beta_{i-1}\over a_{i-1}}y_{i-2}\tag{17} $$ 即 $$ AP_{i-2}w_{i-1}=-{\beta_{i-1}\over a_{i-1}}(r_{i-1}-\gamma_{i-1}p_{i-1}) $$ 代回得 $$ \begin{aligned} \gamma_ip_i&amp;=r_{i-1}-a_{i-1}Ap_{i-1}+{\beta_{i-1}\over a_{i-1}}(r_{i-1}-\gamma_{i-1}p_{i-1})-\beta_{i-1}Ap_{i-1}\ &amp;=\left(1+{\beta_{i-1}\over a_{i-1}}\right)r_i-{\beta_{i-1}\over a_{i-1}}\gamma_{i-1}p_{i-1}\qquad\qquad \qquad\square \end{aligned} $$</p>
<p>最终实现</p>
<div class="sourceCode" id="cb1"><pre class="sourceCode python"><code class="sourceCode python"><a class="sourceLine" id="cb1-1" title="1"><span class="im">import</span> numpy <span class="im">as</span> np</a>
<a class="sourceLine" id="cb1-2" title="2"></a>
<a class="sourceLine" id="cb1-3" title="3"><span class="kw">def</span> CG(A,b,x,n):</a>
<a class="sourceLine" id="cb1-4" title="4">    r <span class="op">=</span> b <span class="op">-</span> A <span class="op">@</span> x</a>
<a class="sourceLine" id="cb1-5" title="5">    <span class="cf">for</span> i <span class="kw">in</span> <span class="bu">range</span>(n):</a>
<a class="sourceLine" id="cb1-6" title="6">        s <span class="op">=</span> r <span class="op">@</span> r</a>
<a class="sourceLine" id="cb1-7" title="7">        <span class="cf">if</span> i<span class="op">==</span><span class="dv">0</span>:</a>
<a class="sourceLine" id="cb1-8" title="8">            p <span class="op">=</span> r</a>
<a class="sourceLine" id="cb1-9" title="9">        <span class="cf">else</span>:</a>
<a class="sourceLine" id="cb1-10" title="10">            p <span class="op">=</span> r <span class="op">+</span> s <span class="op">/</span> _s <span class="op">*</span> _p</a>
<a class="sourceLine" id="cb1-11" title="11">        _s, _p <span class="op">=</span> s, p</a>
<a class="sourceLine" id="cb1-12" title="12">        q <span class="op">=</span> A <span class="op">@</span> p</a>
<a class="sourceLine" id="cb1-13" title="13">        a <span class="op">=</span> s <span class="op">/</span> (q <span class="op">@</span> p)</a>
<a class="sourceLine" id="cb1-14" title="14">        x <span class="op">+=</span> a <span class="op">*</span> p</a>
<a class="sourceLine" id="cb1-15" title="15">        r <span class="op">-=</span> a <span class="op">*</span> q</a>
<a class="sourceLine" id="cb1-16" title="16">        <span class="bu">print</span>(<span class="bu">max</span>(<span class="bu">abs</span>(r)))</a>
<a class="sourceLine" id="cb1-17" title="17">    <span class="cf">return</span> x</a></code></pre></div>
<h2 id="bicg">BiCG</h2>
<p>Let $\mathcal{K}_m=\left&lt;r_0,\cdots,A^{m-1}r_0\right&gt;$, $\mathcal{L}_m=\left&lt;r_0,\cdots,\left(A<sup>{T}\right)</sup>{m-1}r_0\right&gt;$, let $V_m={v_1,\cdots,v_m}$ be a set of basis of $\mathcal{K}_m$, $W_m={ w_1,\cdots,w_m } $ be a set of basis of $\mathcal{L}_m$. We further require that $V_m$ and $W_m$ satisfies $W_m^TV_m=I_m$ . Then $W_m^TAV_m$ and $V_m<sup>TA</sup>TW_m$ are both Hessian, thus $T_m=W_m^TAV_m$ is tridiagonal, $$ T_m= \begin{bmatrix} \alpha_1 &amp; \delta_2 \ \beta_2 &amp; \alpha_2 &amp; \delta_3 \ &amp; \beta_3 &amp; \alpha_3 &amp; \ddots \ &amp; &amp; \ddots &amp; \ddots &amp; \delta_m \ &amp; &amp; &amp; \beta_m &amp; \alpha_m \end{bmatrix}\tag{18} $$</p>
<p>Then we calculate the LU decompostion of $T_m=L_mU_m$, with the diagonal of L<sub>m</sub> equals to 1 $$ T_m=\begin{bmatrix} 1 \ \lambda_2 &amp; 1 \ &amp; \lambda_3 &amp; 1\ &amp; &amp; \ddots &amp; \ddots \ &amp; &amp; &amp; \lambda_m &amp; 1 \end{bmatrix}\begin{bmatrix} \eta_1 &amp; \delta_2 \ &amp; \eta_2 &amp; \delta_3 \ &amp; &amp; \eta_3 &amp; \ddots \ &amp; &amp; &amp; \ddots &amp; \delta_m\ &amp; &amp; &amp; &amp; \eta_m \end{bmatrix}\tag{19} $$ where $$ \begin{aligned} \eta_k&amp;=\alpha_k-{\beta_k\delta_k\over\eta_{k-1}},\quad\eta_0=\alpha_0\ \lambda_k&amp;={\beta_k\over\eta_{k-1}} \end{aligned}\tag{20} $$ Now consider finding $x_m\in x_0+\mathcal{K}<em>m$, such that $r_m=b-Ax_m$ is orthogonal to $\mathcal{L}<em>m$ . It follows immediately that $r_m\in \mathcal{K}</em>{m+1}$ and $r_m\perp\mathcal{L}<em>m $ , then $r_m=kv</em>{m+1}$. Let $x_m=x_0+V_my_m$, then $W_m<sup>Tr_m=W_m</sup>T(r_0-AV_my_m)=0$, i.e. $y_m=T_m^{-1}(\beta e_0)$ , and $x_m=x_0+V_mU_m<sup>{-1}L_m</sup>{-1}(\beta e_0)$ . Let $P_m=V_mU_m^{-1}$ and $Z_m=L_m^{-1}(\beta e_0)$, then by employ that $$ U_m<sup>{-1}=\begin{bmatrix} U_{m-1}</sup>{-1} &amp; -\delta_m\eta_m<sup>{-1}U_{m-1}</sup>{-1}e</em>{m-1} \ &amp; \eta_m<sup>{-1} \end{bmatrix}\ L_m</sup>{-1}=\begin{bmatrix} L_{m-1}^{-1} \ -\lambda_me_{m-1}<sup>TL_{m-1}</sup>{-1} &amp; 1 \end{bmatrix} $$ we can deduce that $$ P_m=[p_1,\cdots,p_m] \quad\text{where}\quad p_m = -{\delta_m\over\eta_m}\underbrace{V_{m-1}U_{m-1}^{-1}e_{m-1}}<em>{p</em>{m-1}}+{1\over\eta_m}v_m \ Z_m=[z_1,\cdots,z_m]^T\quad\text{where}\quad z_m=-\lambda_mz_{m-1} $$</p>
<p>Thus $x_m=x_{m-1}+p_mz_m$ . Now define $P^*_m=W_mL_m^{-1}$ , then $(P^*<em>m)^TAP_m=I_m$ . Now we can derive the BiCG algorithm using following equation $$ r_m<sup>Tr'<em>n=p'<em>mAp_n=0\qquad \text{if}\quad m\ne n $$ where $$ r_n=r</em>{n-1}+z_nAp_n\ r_n'=r</em>{n-1}'+z'<em>nA^Tp'<em>n\ p_n=r</em>{n-1} + \beta_np</em>{n-1}\ p_n'=r_{n-1}'+\beta_n'p'<em>{n-1} $$ then $$ z_n=z_n'={r</em>{n-1}</sup>Tr</em>{n-1}'\over p_{n-1}'^TAp_{n-1}}\ \beta_n=\beta'<em>n={r'^T</em>{n-1}r_{n-1}\over r'^T_{n-2}r_{n-2} } $$ resulting code</p>
<div class="sourceCode" id="cb2"><pre class="sourceCode python"><code class="sourceCode python"><a class="sourceLine" id="cb2-1" title="1"><span class="kw">def</span> BCG(A,b,x,n):</a>
<a class="sourceLine" id="cb2-2" title="2">    r1 <span class="op">=</span> b <span class="op">-</span> A <span class="op">@</span> x</a>
<a class="sourceLine" id="cb2-3" title="3">    r2 <span class="op">=</span> np.array(r1)</a>
<a class="sourceLine" id="cb2-4" title="4">    <span class="cf">for</span> i <span class="kw">in</span> <span class="bu">range</span>(n):</a>
<a class="sourceLine" id="cb2-5" title="5">        s <span class="op">=</span> r1 <span class="op">@</span> r2</a>
<a class="sourceLine" id="cb2-6" title="6">        <span class="cf">if</span> i<span class="op">==</span><span class="dv">0</span>:</a>
<a class="sourceLine" id="cb2-7" title="7">            p1 <span class="op">=</span> r1</a>
<a class="sourceLine" id="cb2-8" title="8">            p2 <span class="op">=</span> r2</a>
<a class="sourceLine" id="cb2-9" title="9">        <span class="cf">else</span>:</a>
<a class="sourceLine" id="cb2-10" title="10">            p1 <span class="op">=</span> r1 <span class="op">+</span> s <span class="op">/</span> _s <span class="op">*</span> _p1</a>
<a class="sourceLine" id="cb2-11" title="11">            p2 <span class="op">=</span> r2 <span class="op">+</span> s <span class="op">/</span> _s <span class="op">*</span> _p2</a>
<a class="sourceLine" id="cb2-12" title="12">        _s, _p1, _p2 <span class="op">=</span> s, p1, p2</a>
<a class="sourceLine" id="cb2-13" title="13">        q1 <span class="op">=</span> A   <span class="op">@</span> p1</a>
<a class="sourceLine" id="cb2-14" title="14">        q2 <span class="op">=</span> A.T <span class="op">@</span> p2</a>
<a class="sourceLine" id="cb2-15" title="15">        a <span class="op">=</span> s <span class="op">/</span> (p2 <span class="op">@</span> q1)</a>
<a class="sourceLine" id="cb2-16" title="16">        x  <span class="op">=</span> x <span class="op">+</span> a <span class="op">*</span> p1</a>
<a class="sourceLine" id="cb2-17" title="17">        r1 <span class="op">=</span> r1 <span class="op">-</span> a <span class="op">*</span> q1</a>
<a class="sourceLine" id="cb2-18" title="18">        r2 <span class="op">=</span> r2 <span class="op">-</span> a <span class="op">*</span> q2</a>
<a class="sourceLine" id="cb2-19" title="19">        <span class="bu">print</span>(<span class="bu">max</span>(<span class="bu">abs</span>(r1)))</a>
<a class="sourceLine" id="cb2-20" title="20">    <span class="cf">return</span> <span class="bu">map</span>(np.array, (X,P1,P2,R1,R2))</a></code></pre></div>
